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ABSTRACT 
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asymptotic expansions that are obtained from integrals with a large parameter. The 
coefficients are defined from recursive schemes obtained from integration by parts. An 
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1. Introduction 

When constructing uniform asymptotic expansions of solutions of differential equations 
or of functions defined by integrals, usually a difficulty arises when the coefficients of the 
expansion are constructed. As shown in Olver (1974) for the Airy-type expansions of 
Bessel functions, recursion relations for the coefficients can be obtained for the case that 
the expansion is obtained by using a linear second order differential equation. 

In many publications this method has been used, for example in Olver (1959) and 
Dunster (1989), and for expansions involving Bessel functions or parabolic cylinder func- 
tions similar results are available. Having such a recursion relation for the coefficients does 
not always give the possibility to obtain analytic expressions of a number of coefficients, 
because the recursion involves integrals of previous coefficients together with a function 
that is not easy to handle. Sometimes the coefficients can be explicitly expressed in terms 
of coefficients of simpler expansions because different types of expansions may be valid 
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in overlapping domains. See for the Bessel functions the relations in Olver (1974), page 
425, Exercise 10.3 or Abramowitz & Stegun (1964), page 368, formula 9.3.40. 

For special functions usually the same type of uniform expansions can be obtained 
by using integral representations of the functions. Sometimes, in a particular problem, 
the integral is the only tool available for constructing uniform expansions. By using 
transformations of variables in the integrals, these representations can be transformed 
into standard forms for which an integration by parts procedure can be used to obtain 
expansions in terms of, for example, Airy functions. 

Although it is usually not possible to derive recursion relations for the coefficients 
obtained in this way, in all cases for special functions known so far, it is possible to 
construct a number of coefficients, and only because of the complexity of the problem, 
which implies limitations with respect to available computer memory when doing symbolic 
computations, there is an upper bound for this number. An advantage of the differential 
equation approach is the possibility to construct realistic and sharp error bounds for the 
remainders in the expansions; similar bounds cannot be obtained in the approach based 
on integral representations. 

In this paper we use integral representations and give Maple algorithms for construct- 
ing the coefficients in uniform asymptotic expansions involving Airy functions. First we 
describe how to obtain the coefficients for a general case. For an application we obtain the 
coefficients for the case of a special function called parabolic cylinder function. Straightfor- 
ward computations are often complicated by appearance of algebraic roots in the output 
or intermediate expressions. These algebraic roots can be avoided by replacing some pa- 
rameters with algeraic expressions in suitable new variables. In the example of parabolic 
cylinder function we avoid computations with algebraic roots by using variable u (instead 
of t), and at the end we simplify the output by introducing variable £. In the last section 
we give the Maple code used for this example. 

2. Airy-type asymptotic expansions 

We consider integrals of the form 



where the contour is starts at infinity with phi = — it/ 3 and returns to infinity with 
phi = 7r/3. We assume that the function f(t) is analytic in the neighbourhood of the 
contour. The parameter z is large positive number, and rj is also assumed to be real. 
Extension to complex values of the parameters is possible, but this will not be discussed 
in this paper. 

In the case f(t) = 1 we obtain the Airy function (Temme (1996), page 101) 




(1) 




(2) 
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For more general functions / the asymptotic expansion of F v (z) can be given in terms 
of this Airy function. The asymptotic feature of this type of integral is that the phase 
function (j)(t) = |t 3 — rjt has two saddle points at ±^/r] that coalesce when r/ — > 0, and it 
is not possible to describe the asymptotic behaviour of F v (z) in terms of simple functions 
when rj is small. When the parameter rj is positive and bounded away from 0, one can 
perform a saddle point analysis on (1) and use a conformal mapping (f>(t) — (f)(-y/rj) = \u 2 
with the condition u(^/rj) = 0. We obtain 

-i />ioo 

F v (z) = — / e$™*g(u)du, 
2m J_ iQO 

where g(u) = f(t) dt/du, with dt/du = u/(t 2 — rj), which is regular at the positive saddle 
point, but not at the negative saddle point. It follows that, when ry becomes small, a 
singularity due to dt/du in the it— plane approaches the origin, and an expansion of dt/du 
at u = will have coefficients that become infinite as rj — ► 0. Hence, by using the standard 
saddle point method we obtain an expansion that is not uniformly valid as rj — ► 0. 

A modification of the saddle point method is possible by taking into account both 
saddle points. We give an integration by parts procedure that is a variant of Bleistein's 
method introduced in Bleistein (1966) (for a different class of integrals), and that gives 
the requested uniform expansion. 

We assume that / is an analytic function in a certain domain G and write 

f(t) = a + Pot+(t 2 -vMt), (3) 



where 

ao = l[f(Vv) + f(-Vv)], fo = 2j=[HVv)-H-V=v)]- (4) 

Clearly ao — ► /(0), (3q — > /'(0) as rj — > and the following Cauchy integral representations 
hold 

"0 = 7T- / 1 A) = 7T~ / 1 ds > 

2m J q s Tj 2m J c s z - rj 

where the contours of integration C encircle the points t and/or ±i/rj. Upon substituting 
(3) in (1), we obtain 

F v {z)=z-*M (pz*}ao-z-*M' (yz*} (3 + ^- J e ^¥ 3 -^){t 2 - rf)g{t)dt. 

An integration by parts gives 

F v (z) = z-^Ai (pzi) a - z~*M' (pz*} O - ^- J e z ^-^h{t) dt, 
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where fi(t) = g'(t). Repeating this procedure we obtain the compound expansion 

OO OO q 

F v (z) ~ z-*Ai (^*) £(-1)"^ " ^ §Ai ' § ) E(-!) n ^' ( 5 ) 



n=0 n=0 



where the coefficients a n , /3„ are defined as in (4) with the function / replaced with f n , 
which in turn is defined by the scheme 

f n+1 (t) = g' n (t), f n (t) =a n + (3 n t + (t 2 - r))g n (t), (6) 

with n = 0, 1, 2, . . . and fo(t) = f(t). The expansion in (5) is valid for large values of z 
and holds uniformly with respect to r] in a neighbourhood of the origin. A more precise 
formulation can be given, but more information can be found in the literature; see Olver 
(1974) and Wong (1989). 

The functions f n (t) defined in (6) can be represented in the form of Cauchy-type 
integrals. We have the following theorem. 

Theorem 1. Let the rational functions R n (s,t,rf) be defined by 

R (s,t,rj) = — R n+1 (s,t,ri) = ^ ^ - ^-R n (s,t,r)), n = 0,1,2,..., (7) 

where s, t, rj G (D, s / t, s 2 / rj. Let f n (t) be defined by the recursive scheme (6), where fo 
is a given analytic function in a domain G. Then we have 

fn(t) = J Rn(s,t,ri)f (s)ds, 
where C is a simple closed contour in G that encircles the points t and ±y/rj. 
Proof. The proof starts with 

fn(t) = J R Q (s,t,r])f n (s)ds, 

and in this representation the recursion relation (6) for the functions f n is used. More 
details can be found in Olde Daalhuis & Temme (1994). g 

For the coefficients a n ,[3 n we have a similar representation: 

an = ^~j A n(s,V)fo(s)ds, Pn = ^~j B n(s , 1]) f (s) ds , (8) 

where C is a simple closed contour in G that encircles the points ±yfrj and where A n (s,t) 
and B n (s, t) follow the same recursion (7) as the rational functions R n (s, t, rj), with initial 
values 

s 1 
A o(s,v) = -o > B o{s-,rf) = — . 

S z — Tj S z — f] 
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We see that the coefficients a n ,(3 n that play a role in the expansion (5) are well 
defined from an analytical point of view. However, from a computational point of view 
it may be quite difficult to evaluate the coefficients. For a simple rational function like 
fo(t) = l/(i + l) the computations are rather straightforward, and we can even use residue 
calculus to evaluate the integrals in (8): 

a n = -A n (-l,rf), p n = -B n (-l,rf). 

The first few values are in this case 

1 i 

«0 = 7' A) 



1 ' — -.i 

i] — 1 r/ — 1 

ttl (77-1)3' Pi (r/-l) 3! 

«2 = -4- 777, /3 2 = 2- 



(77 -1)5' " (,,-1)5' 
2^ + 217/ + 7 _ , n r] + 2 



(r?-l)7 ' ^" fa -1)7' 
OQn 7? 2 + 4?? + l ?7 2 + 1977 + 22 

a4 = " 280 ( V -D* ' ^ = 4 ° („-!)» ■ 

OQn ^ 3 + 29T7 2 + 6577 + 13 on 2r/ 2 + M77 + 11 
«5-280 (^_i)n ' /% = -1120 7— ^ n . 

For a more complicated or general function /q (t) even computer algebra manipulations 
give complicated expressions which are very difficult to evaluate. In the next section 
we develop an algorithm for computing the coefficients a n , (3 n when the values of the 
derivatives of fo(i) at t = iy^ are available. 

3. How to compute the coefficients a n ,f3 n 

We explain how the coefficients a n ,(3 n of (5) can be computed. To avoid the square roots 
in the formulas we replace 77 with b 2 , and we write (6) in the form 

fo(t) = f(t), f n+1 (t) = g' n (t), f n (t) =a n + p n t + (t 2 - b 2 )g n (t), 

for n = 0, 1, 2, . . .. We assume that the function / is analytic in a domain G, that the 
series expansions used in this section are convergent in G, and that the points ±6 are 
inside G. Furthermore, we assume the coefficients p^\p^ of the expansions 

00 00 

/(*) = Epi 1) (*-6) fc . n-t) = Y,p i k ) ( t - b ) k ( 9 ) 

k=0 k=0 



are available. 
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Theorem 2. Algorithm. 

Let coefficients f%, f% be defined by 



J k r> 



(1) i (2) 
Pk +Pk 



f° = - 
> Jk 2 



(1) (2) 
Pk ~Pk 



k = 0,1,2... 



and coefRcents f^ ,e by the recursion 

bf° k ' e = fk ~ fk'-v k>0, 
with f^l = 0. Next, define coefficients 7^,5^ by 

7o = fo, <5o = fo' e , 



and for k > 1: 



7fc = J^ 



■l) fc -^-(2fc-j-l)! 



(2b) 2k ~3 k\(k- j)\ 



f e 
j 1 j 



A (-i)fc-jj(2fc-j-l)! 



^ (26) 2 *-J fc!(fc- j)! 
Finally, let for n > coefficients , be defined by the recursion 



(n+l) 



(2fc + 1)4^ + 2^ + 1)^2, 



(n) 



2(fc+ 1)7^2, fc = 0,1,2,... 



(10) 



(11) 



with 7^ = 7fc,5^ = 5k- Then the coefficients os n ,(3 n of expansion (5) are given by 



a, 



= 7?°, Pn = ^\ n>0. 



Proof. The coefficients /| , f% occur in the expansions 

oo oo 
fe(t)=Y.ft{t-b) k f (t) = ^fk(t-b) k , 



k=0 



k=0 



where f e (t), f Q (t) are the even and odd parts of /: 



fe(t) = l[f(t) + /(-*)], /o=^[/(*) "/(-*)], 



and the coefficients /^' e occur in the expansion 

1 oo 

-jo{t) = Y J f°u e {t-b)\ 



k=0 
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The coefficients jk^k occur in the expansion 

oo oo 

f(t) = E ^(t 2 - b 2 ) k +^M« 2 - b 2 ) k . (12) 

fc=0 k=0 

Observe that 

oo oo 

f e (t) = E 7fc(i 2 - ft 2 )", /<>(*) = tE^ 2 " 62 ) fc ' 

fc=0 fc=0 

and we will verify the first relation of (10). We write 

dz 



z k+r 

where the contour is a small circle around the origin. Also, 

2tdt 



Ik 



5S /'•<*> 



(t + b) k + 1 {t-b) k + 1 ' 

where the contour is a small circle around t = b. 

Substitute the expansion f e (t) = YlJLo ~~ &) J '- Then, 

k 



Expand 



_ 1 f 2tdt 

lk ~^J j 2^1 J (t + b) k + 1 (t-b) k + 1 -i' ( } 



9 , oo 

- T =J2<lrn(t-br. (14) 



(t + b) k + 

' m=0 

We find, by using (Temme (1996), p. 108) 

n=0 ' n=0 V y 

_ / ^ m (k-m)(k + m-l)\ 
Qm [ } {2b) k + m m\k\ ' 

When we use (14) in (13), we only need q m with m = k — j. This gives the first result of 
(10). The proof for 5k is the same, because (l/t)f (t) is again even. 
The coefficients 7^ , 5^ are used in 

00 00 

/»(*) = E - ^) fc + * E *£V - 62 ) fc > 

fc=0 fc=0 

and the recursions in (11) are easily verified, as is the final relation 

«n = 7o ) Pn = 70 J n > °- 
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The first few values are of the coefficients 7fc,<5fc of the expansion in (12) are 

7o = /o) <*o = ^/o > 

7i = ^/r, ^1 = ^3^/1-/0), 

72 = 2^Wf " /f), 5 2 = ^(26 2 / 2 ° - 36/f + 3/ °), 

and we observe, as in (10), negative powers of b. From a computational point of view, 
this may cause numerical instabilities, because the coefficients are analytic functions of b 
at b = 0. For example, taking f(t) = l/(i + 1) again, we obtain 

7fc = (1 _ fo2)fc+i' 5fc = ~(i _ 62)fc+i' & = 0,1,2,..., 
which follows from 

1 l-t l-t ^ (t 2 -fe 2 ) fc ^ (t 2 -fe 2 ) fc 



t + l l-t 2 (1 - b 2 ) - (t 2 - b 2 ) ^ {l-b 2 )^ 1 ^ (1 - &2)fc+i • 



From the representations in, for example, (10), we conclude that if we apply the al- 
gorithm for computing the coefficients a n ,(3 n of expansion (5), starting with numerical 
values of the coefficients p^\p^ of (9), we may encounter numerical instabilities when 
b is small. For this reason, it is important to use exact values of p^\p^\ and computer 
algebra is of great help here. In the next section we consider a non-trivial case in which 
obtaining the exact values of the coefficients p^\p^ of (9) also needs special care. 



Remark 1. In order to compute the coefficients a n ,/3 n for n = 0,1,..., N from the 
relation 

and the recursion in (11), we need the starting values for this recursion 7^,(5^ for k = 
0, 1, . . . , 2N. Hence, as follows from (10), we also need p^\p^\ k = 0, 1, ... , 2N in the 
expansions in (9). 



4. Application to parabolic cylinder functions 

Weber parabolic cylinder functions are solutions of the differential equation 

g-(I* 2 + a) y = 0. (15) 

Airy-type expansions for the solutions of this equation can be found in Olver (1959), and 
are obtained by using the differential equation. In this section we show how to obtain an 
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integral representation like (1), and how to apply the algorithm of the previous section for 
deriving an Airy-type asymptotic expansion. 

A standard solution of (15) is the integral (see formula (19.5.4), page 688 in A&S) 

U{a,x) = ^= / e -*«+*»V«-i ds, (16) 
W2ir Jc 

where the contour C is a vertical line in the complex plane with 5fts > 0. 
We consider large negative values of a, and use Olver's notation 

a = —TjfJ 2 , x = fitV2. 

Changing the variable of integration by writing s — > [is/ y/2, we obtain 



where 



(s) = \s 2 - 2st + lns, z = \\x 2 . 



The saddle points are obtained from the equation <j)'(s) = 0, that is, from 

s 2 -2st + l 
s 

which gives two solutions 



t±Vt 2 - 1. 



The saddle points coalesce when t — > ±1. Observe that in the new variables the differential 
equation (15) transforms into 

which has turning points at t = ±1. 

A transformation into the standard form (1) can be obtained by writing 

(f)(s) = lw 3 -rjw + A, (17) 

where rj and A have to be determined and do not depend on t. A transformation into 
the cubic polynomial is first considered in Chester et al. (1957). For further details on 
the theory of this method we refer to Olver (1974), Wong (1989), and Olde D. & T 
(1994). 

The parameters r\ and A are obtained by assuming that the saddle points s± in the 
s— variable should correspond with the saddle points w± = ±^/r\ in the w— variable. We 
write 

t = cosh6*, which gives s± = e ±6 , (18) 



10 



assuming for the time being that 9 > 0. We obtain the equations 

i e+2 , _ 2e+e cQsh9 + e = _ 2^3/2 + ^ 

\ e ~™ _ 2e - e coshtf - 9 = +f rfl 2 + A, 

from which we derive 

f?? 3 / 2 = sinh2#-2#, A = -\ - cosh 2 9 = -\ - t 2 . (19) 

By using these values of rj and A the w— solutiuon of the equation in (17) is uniquely 
defined. Namely we use that branch (of the three solutions) that is real for all positive 
values of s, and s > correponds with w £ IR. 

After these preparation we obtain the standard form (cf. (1)) 

U (-±m 2 , »tV2) e~* A = y^e^ 2 ' 2 (f ) ^ F v (z), (20) 

where 

w = ^// <i "*-''" , /(w)<;w, /w = ^- 

Taking into account the mapping in (17), we have 

ds w 2 -b 2 .. , ^ w 2 - b 2 j2 

d»- V-a.+i - /w = ^-2 t » + r < 21 > 

As explained in the previous section, for the computation of the coefficients a n ,(3 n , we 
need the coefficients p^\p^ of the expansions (cf. (9)) 

oo oo 

f{w)=^p^(w-b) k , f(-w)=Y J P i *\™-b) k . (22) 

k=0 k=0 

It turns out that p^ = p^ . Indeed, consider the expansions: 

oo oo 

s = s+ + J2s+(w-b) k , s = s _ + J2 Sk -(w + b) k . (23) 

k=l k=l 

Using the expression of ds/dw in (21) and l'Hospital rule we obtain 



2b , , \/bs+ I bs. 

s_l_— — so that 



'2(s+-t)sJ 1 (t 2 -l)^ V sinhtf 

The square root has the plus sign because ds/dw is positive if w € IR, as follows from the 
first relation in (21) and the properties of the mapping. From the expression (24) for f(w) 
we obtain: 



P { o> = f(b) 



y/s^ V sinh# 
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Analogously, 



bs. 



sinhe and V » ) = f{ - h) = ^e 



- J 1 ) 



PV- (25) 



In order to avoid expressions with algebraic roots in the computations, it is convenient 
to consider expansions like (22) for the function f(w) = f(w)/f(b). The corresponding 
coefficients we denote by p^ and p^K Besides, to avoid algebraic roots in the expansions 
of (23) we replace t by a new variable 

r-,'t-l\* , 46 2 +n 4 

u = V2o , so that t 



t + lj ' 46 2 -u 4 ' 
Then 

_ 2b + u 2 + _2b + u 2 
S+ ~ 2b^' Sl ~ ~^u~' 

Other coefficients s\ can be obtained by deriving a reccurence relation for them from 
the differential equation in (21). They are rational functions in u and b. The coefficients 
can be obtained from the corresonding s~l by changing the sign of both u and b. In 
particular, 

_ 2b -u 2 _ _ 2b -v 2 
S ~ ~ 26+^?' Sl ~ ~^hT ' 

Further, the coefficients p^ and p^ can be computed using 

\J s aw aw 

Recall that y/s satisfies the differential equation 2sdS/dw = Sds/dw. It is convenient to 
compute the power series (in w — b) solution S+ of this equation with S+(b) = Au/(2b—u 2 ). 
Then f(w) = dS+/dw and the coefficients p^ are obtained easily. The coefficients p^ 
can be obtained by changing the sign of both b and u in the expression for {—\) k p^\ 

Application of the algorithm of the previous section gives the coefficients a>j, (3j for 
the expansion of f(w), and these coefficients are rational functions in b and u. We write 
them in a more compact form as rational functions in rj = b 2 and 

u 4 + 4b 2 bt 



s Au 2 " v/P^T' 
Then first few coefficients in the expansion (5) are: 

a = 1, Po = 0, 

1 5£ 3 - 6r/£ - 5 
01 " 48' A " 48^ ■ 
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385£ 6 - 924r/£ 4 + Q84t] 2 ( 2 - 143r/ 3 + 70£ 3 - 84??£ - 455 



a 2 



4608 f] 



3 



03 



a Pi « 2 2021 

& = 48' a3 " 48" 34560 ai ' 
425425£ 9 - 15315307?^ + 2040012r/ 2 £ 5 - 28875£ 6 - 1189005r/ 3 £ 3 + 69300^ 4 

3317760 ry 5 

259110?7 4 C - 51300r/ 2 £ 2 + 28875£ 3 + 10725r/ 3 - 34650^ - 425425 



3317760 r? 5 



B - A 2021 fl 
^ 4 " 48 " 34560^- 



The linear relations between the coefficients follow from expansion (8.11) in Olver(1959), 
where both power series factors of Ai and Ai' contain only even powers of our z (in Olver's 
notation, z = \^ 2 )-, but the whole expansion is multiplied by function g(z) with known 
asymptotics. Olver also notes that coefficients in the asymptotic expansion of U(a, x) 
in terms of Airy functions can be lineraly determined from the asymptotic expansion 
(of the same function) in terms of elementary functions; see formulas (8.12), (8.13) in 
Olver(1959). 

The coefficients a n , (3 n are analytic functions at 77 = and we can expand them in 
Maclaurin series. The first few coefficients are expanded as follows: 

9 7 1359 2 7 3 _ 152723 4 3997 5 

f ~560 + 1800 71 ~ 1078000^ + 16250^ ~ 1018710000 ^ + 75968750 V + '"' 
199 6849 737 2 46711 3 975823 4 

" 2 " ~ 115200 + 4928000 V ~ 1040000^ + 142560000^ ~ 6806800000^ + "" 

The radius of convergence equals (37r/2) 2 / 3 = 2.81.... This number follows from the 
singularity of the mapping given in (19), with 6 defined in (18). The mapping is singular 
at t = -1. 



5. Maple code 

General case. For an input one has to (re)define functions AiryPw and AiryPm specifying 
the coefficients in (9). Output is given by functions AiryAlpha and AiryBeta, which return 
the coefficients in (5). For convenience, one may rename the global variable AiryB using 
alias. 

AiryAlpha:= proc( n ) normal(AiryGamma(n,0)) end: 
AiryBeta:= proc(n) normal(AiryDelta(n,0)) end: 
AiryGamma:= proc( n, k ) 
if n=0 then AiryC(k) 

else factor( (2*k+l)*AiryDelta(n-l,k+l)+2*AiryB' 1 2*(k+l)*AiryDelta(n-l,k+2) ) 
fi 

end: 

AiryDelta:= proc(n, k ) 
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if n=0 then AiryD(k) 

else 2*(k+l)*AiryGamma(n-l,k+2) 

fi 

end: 

AiryC:= proc( k ) local j; 
if k=0 then AiryFe(O) 
else factor( 

sum( , H)' 1 (H)*j/(2*k-j)*binomial(2*k-j 1 k)/(2*AiryB)' 1 (2*k-j)*AiryFe(j)' 1 'j'=l..k ) ) 

fi 

end: 

AiryD:= proc( k ) local j; 
if k=0 then AiryFoe(O) 
else factor( 

sum( , H)' 1 (k-j)*j/(2*k-j)*binomial(2*k-j,k)/(2*AiryB)' 1 (2*k-j)*AiryFoe(j)', 'j'=l..k ) ) 

fi 

end: 

AiryFoe:= proc( k ) 
if k<0 then 

else expand( (AiryFo(k)-AiryFoe(k-l))/AiryB ); 
fi 

end: 

AiryFe:= proc(k) (AiryPw(k)+AiryPm(k))/2 end: 
AiryFo:= proc(k) (AiryPw(k)-AiryPm(k))/2 end: 

To find the coefficients of the expansion of parabolic cylinder function U (a, x) (up to 
the multiple in (25)) one has to assign 

AiryPw: = ParCyPw; AiryPm:= ParCyPm; 
The global variables are AiryB, ParCyU, ParCyXi, they correspond to variables b,u,£ in 
the text. The coefficients in b and u would be returned by AiryAlpha and AiryBeta, and 
coefficients in b and £ — by ParCyAlpha and ParCyBeta. 

alias( ParCyUa=RootOf( z"4-4*ParCyXi*z"2+4*AiryB"2, z)): 
# Algebraic relation between ParCyU and ParCyXi 

ParCyAlpha:= proc(k) factor( evala(subs(ParCyU=ParCylla,AiryAlpha(k))) ) end: 
ParCyBeta:= proc(k) factor( evala(subs(ParCyU=ParCyUa,AiryBeta(k))) ) end: 
ParCySw:= proc( k ) option remember; local T, s, a, w; 
if k=0 then (2*AiryB+ParCylT2)/(2*AiryB-ParCyir2) 
elif k=l then (2*AiryB+ParCylT2)/2/ParCyU 
else s:= sum('a[i]*wT, 7=0. .k); 
T:= coeff( expand( (s"2+2*(ParCyU"4+4*AiryB"2)/(ParCyU' 1 4-4*AiryB"2)*s+l) 

*diff(s,w)-w*(w+2*AiryB)*s ), w, k); 
sort( factor( solve( subs( seq(a[i]=ParCySw(i),i=0..k-l), T), a[k])), ParCyU) 

fi 

end: 

ParCySm:= proc(k) subs( ParCyU=-ParCyU, AiryB=-AiryB, ParCySw(k) ) end: 
ParCySqrtS:= proc( k ) option remember; local j; 
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if k=0 then 4*ParCyU/(2*AiryB-ParCyU"2) 

else factor( sum('(3/2*j-k)*ParCySw(j)*ParCySqrtS(k-j)', 'j'=l..k) /ParCySw(0)/k ) 
fi 

end: 

ParCyPw: = proc(k) (k+l)*ParCySqrtS(k+l) end: 

ParCyPm:= proc(k) (-l)"k*subs(AiryB=-AiryB, ParCyU=-ParCyU, ParCyPw(k)) end: 
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